Electrostatics and aggregation: how charge can turn a crystal into a gel 
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The crystallization of proteins or colloids is often hindered by the appearance of 
aggregates of low fractal dimension called gels. Here we study the effect of elec- 
trostatics upon crystal and gel formation using an analytic model of hard spheres 
bearing point charges and short range attractive interactions. We find that the chief 
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electrostatic free energy cost of forming assemblies comes from the entropic loss of 
counterions that render assemblies charge-neutral. Because there exists more accessi- 



ble volume for these counterions around an open gel than a dense crystal, there exists 
an electrostatic entropic driving force favoring the gel over the crystal. This driving 
force increases with increasing sphere charge, but can be counteracted by increasing 

co 

>• counterion concentration. We show that these effects cannot be fully captured by 

in 

pairwise-additive macroion interactions of the kind often used in simulations, and we 
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show where on the phase diagram to go in order to suppress gel formation. 
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I. INTRODUCTION 

Crystallizing proteins or colloids is difficult in part because the growth of kinetic ag- 
gregates called gels can compete with crystal nucleation 1-3 . Considerable insight into the 
competition between crystallization and gelation has been obtained from models of spheri- 
cal particles in implicit solvent with pairwise-additive interactions, modeling both the forces 
(hydrogen bonds, hydrophobic effects etc.) that drive association, and the electrostatic re- 
pulsions 4 that act to suppress association 5-9 . Such models have shown how the phase diagram 
varies with interaction strength and range, and have revealed the presence of kinetically- 
arrested gel- and glass states 10 . While pairwise-additive 'screened' potentials capture much 
of the important physics of electrostatics, they neglect explicit counterion degrees of free- 
dom. Here we show that accounting for counterion entropies reveals that electrostatics does 
more than simply prevent association: it can induce a thermodynamic driving force that 
favors a gel over a crystal. 



II. MODEL 

We consider a collection of hard spheres (macroions) of radius a in aqueous solution with 
added salt (we set a = 1.6 nm, appropriate for a small protein). Spheres carry a point 
charge q at their centers, and interact with square-well attractions of depth — e and range 
2a(A — 1). We are interested in three phases: the solution phase of dissociated spheres, 
a close-packed three-dimensional crystal, and a kinetically-stabilized gel that we model as 
a one-dimensional rod of associated spheres. We describe a crystallization window that is 
bounded on one side by the thermodynamic stability of the crystal, and on the other by 
the appearance of the kinetically favored gel. We approximate the kinetic threshold for 
gelation as the regime in which the chemical potential of a linear aggregate drops below 
that of the solution; since one-dimensional objects lack a nucleation barrier, we expect 
gelation to preempt crystallization in this regime. This threshold reproduces qualitatively 
the gelation boundaries determined from numerical and experimental studies of colloids and 
proteins 3,7 ' 11-13 while allowing efficient computation of electrostatics. 

We write the chemical potentials of the zero-, one- and three-dimensional phases (shown 
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FIG. 1. Counterion entropies dominate the electrostatic free energies of d-dimensional assem- 
blies. We plot per-particle electrostatic energies, Ud, entropies, Sd, and free energies fd for each 
d-dimensional assembly (pictured left; the case d = 2 is included for comparison only). We have 
assumed a macroion charge q = 5e and have considered four different salt concentrations c. The 
electrostatic free energy cost for establishing each assembly is dominated by the entropy loss upon 
counterion confiment, and not by macroion-macroion repulsion energies. 



in Fig. 1) as 



jiq = lnp + p(4 — 3p)(l — p) 2 solution; 

p 1 = -e-A;BTln(A 3 -l) + A/r s gel; 
p 3 = -ze/2 - 3k B T ln(A - 1) + A/ 3 0S crystal. 



(1) 



OS 

d 



These expressions account for free energies of binding and electrostatic free energies A/ ( 
associated with each d— dimensional assembly {d — refers to solution-phase monomers, 
d = 1 refers to chain-like aggregates and d = 3 describes the crystal phase). The chemical 
potential po for the solution phase contains the Carnahan-Starling expression for the free 
energy of hard spheres as a function of p, the volume fraction of spheres in solution 14 . 
The expression \i\ for the rod phase accounts for the fact that each particle makes two 
pairwise contacts of energy — e. It also contains an entropic term that accounts for the 
fact that a new particle can add to the end of a rod anywhere within a shell of volume 
~ (2aA) 3 — (2a) 3 around the existing end particle. This contribution accounts for both 
bending and vibrational entropies (we neglect unimportant constant terms). We also neglect 
rod branching 15 , which at the level of electrostatics amounts to assuming the characteristic 
linear extent of a gel is long compared to the Debye length. The chemical potential p 3 of 
the close-packed crystal accounts for the energetic contacts between each particle and its z 
nearest neighbors (we set z = 12), and the fact that each particle can vibrate a distance of 



order 2a(A — 1) in each direction 6 . We set A = 1.2 in our calculations, motivated by the 
recognition that proteins typically possess interactions of short range 6,16-18 . 

To these expressions we add the electrostatic free energies A/J s of each rf-dimensional 
assembly. We calculated these free energies using the Poisson-Boltzmann (PB) equation 

V 2 $ d = sinh($ d ), (2) 

which we have written in terms of scaled spatial coordinates x = nr and a dimensionless 
potential $d = eipd/k^T. k~ x is the Debye length, given by k 2 = 2e 2 c/ (ek-^T), where c is 
the concentration of ions in a reservoir of solvent in osmotic equilibrium with the system. 
We assume our charged particles to be in aqueous solution with positive and negative salt 
ions of local concentration c d (r) = ce =F * d( - r - ) . To determine the potential $^ we solved 
Eq. (2) for given assembly geometry. To simplify the solution of the PB equation we 
approximated the rod as a smooth cylinder. For the crystal we used a cell model 16 ' 19 ' 20 . We 
have also considered a two-dimensional assembly to gain insight into the scaling behavior of 
the electrostatics, although we stress that one would need to consider an anisotropic particle- 
particle interaction to stabilize such a structure (see e.g. 21,22 ). The electrostatic calculations 
are presented in detail in the appendix. 

We then computed the electrostatic energy Ud and entropy S d associated with each as- 
sembly. The energy Ud stored in the electric field is 
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U d =-J tfV(Vr^) , (3) 

where e = 80e is the permittivity of water (e is that of free space). The entropy S d 
quantifies the cost of confining counterions and excluding coions from the screening layer 
associated with a rf-dimensional assembly 23 . The entropy cost per unit volume of maintaining 
a screening layer of ion concentration c s in thermal contact with a reservoir of ions at 
concentration c is 

-S/ks = fdc' In (c'/c) 

= c s ln(cg/c) -c s + c. (4) 

The electrostatic free energy can then be written F d = Ud — TSd = f d 3 rJ : d, where the free 
energy density is 

^ = |(V r ^) 2 + A;BT{c+ln(c+/c) 

+^hi(c d /c)-c+-c d +2c}. (5) 



Note that the minimization of Eq. 5 with respect to ionic concentration results in Eq. 2 24 . 
Finally, to compute the electrostatic free energy cost of assembling spheres into each 
d— dimensional structure, we computed the free energy per particle in the assembly, fd = 
Fd/N, and subtracted from this the electrostatic free energy of the monomeric building 
blocks in solution (determined by solving the PB equation for an isolated sphere). The 
result is A /J 3 . 

III. RESULTS 

A. Counterion entropy dominates the electrostatic free energy of assembly. 

Fig. 1 shows the electrostatic free energy fd of a rf-dimensional assembly, and its compo- 
nents, the electrostatic energy Ud and entropy Sd (the case d — 2 is included for comparison, 
but it will not appear in our phase diagram). We see that the energy per particle is non- 
monotonic in aggregate dimensionality: the energy change upon condensation from the free 
particle state to the crystalline state is negative, because the macroion-macroion repulsion 
is overcome by the attraction between macroions and the counterions confined within the 
crystal's cavities (assumed to be small compared to the Debye length). By contrast, the 
electrostatic entropy cost increases monotonically with aggregate dimensionality. Further- 
more, these entropies are much larger in magnitude than their energetic counterparts, and so 
the electrostatic free energy cost becomes larger as the structure increases in dimensionality. 
Although it seems intuitively obvious that assembling charged spheres into a crystal should 
meet with an electrostatic free energy cost, it is not obvious, a priori, that the origin of this 
cost is entropic, and not energetic 25 . 

Fig. 1 reveals that the electrostatic entropy alone is a reasonable approximation of the 
electrostatic free energy cost of assembly, a fact that has been long appreciated in polymeric 
systems 26 ' 27 . The behavior of the entropy with d can be understood from the following simple 
argument. Electro neutrality requires that q counterions be confined near each macro ion. 
Counterions are free to extend a distance k~ 1 in each of the 3 — d dimensions extending away 
from the assembly, but are confined within a distance of order the particle size a in each of 
the d dimensions of the assembly. Thus the volume v localize within which q counterions must 
be localized scales as v localize ~ a d n d ~ 3 , and so the entropic cost per macroion of building the 



screening layer scales as —Sd/k ~ qln(c(an) d /qK 3 ), or 

—Sd/k ~ qd\n(aK). (6) 

This cost increases with d. 

B. Electrostatic interactions are not pairwise additive. 

One immediate consequence of the fact that counterion entropies dominate the free energy 
cost for assembling charged particles is that effective macroion interactions are not in general 
pairwise-additive: the free energy difference A/J 3 depends on the volume of space accessible 
to counterions, but does not necessarily vary linearly with the number of contacts made by 
particles in each assembly We can define A31 = zA/i s /2A/| s as a measure of the relative 
electrostatic cost of making a macroion-macroion contact in a gel versus a crystal. When 
A31 is unity there is no difference in this cost, the electrostatic effects can be adsorbed into 
an effective pairwise interaction, such as a Yukawa potential 28 . However, when A 31 7^ 1 the 
electrostatic free energy is non- additive, rendered so by counterion degrees of freedom. 

In Fig. 2 we plot A31 as a function of salt concentration c and particle charge q. A31 is in 
general less than unity, indicating that it is easier, per macroion contact, to fit counterions 
around a gel than within a crystal. A31 is also non-monotonic: in principle, at large salt 
concentration strong electrostatic screening will ensure that macroion sites are independent, 
with electrostatic free energies becoming pairwise-additive (as seen in the plot, though, such 
concentrations can be unattainably high). At low salt concentration A 31 again approaches 
unity. This can be understood by considering Eq. 6 for a cubic packing arrangement. In 
this case d is proportional to the number of nearest neighbors and each additional pair of 
neighbors reduces the volume accessible to the counterions by a factor na. This results in 
an electrostatic free energy linear in the number of nearest neighbors. However, at finite 
salt concentration the volume occluded by the macroion (neglected in Eq. 6) becomes a 
significant perturbation to localize, and the approximation of pairwise additivity fails. 

C. Model phase diagram: how does one avoid gelation? 

We turn now to calculation of the model's phase diagram. We found the solution-crystal 
coexistence line by setting /x = /i 3 , and estimated the nonequilibrium solution-gel coexis- 
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FIG. 2. Relative electrostatic cost of making a macroion-macroion contact in a gel versus a crystal, 
A31, as a function of salt concentration c and macroion charge g. Only when A31 is unity is the 
electrostatics of assembly effectively pairwise-additive, which in general is not the case. 



tence line by setting \Iq — \X\. This line does not describe an equilibrium coexistence: what 
we mean by it is the location at which there exists a driving force for monomers to aggre- 
gate in a disordered way. The crystal may be thermodynamically stable, but at moderate 
degrees of supercooling or supersaturation it is separated from the solution phase by a free 
energy barrier, and may require considerable time to appear. By contrast, there exists no 
free energy barrier to the formation of a linear aggregate, and so we expect a gel to form 
readily below the nonequilibrium solution-gel coexistence line. 



Our phase diagram in the temperature-density plane is shown in Fig. 3(a), for fixed 
salt concentration and for two macroion charges. The close-packed crystal is stable below 
the gray solution-crystal coexistence lines, but we expect a driving force for gelation below 
the red solution-gel coexistence lines. This behavior is similar to that seen in simulations 
of spheres with short range attractions 12 . The resulting 'window' of crystal stability and 
accessibility (shaded region between the red and blue lines) is narrow, and becomes narrower 
as macroion charge increases. This is a key finding of our study. As noted before, the reason 
for this narrowing is the increasing cost of confining counterions within a close-packed crystal 
relative to around an open gel. 

D. High salt concentrations promote crystallization. 

We have seen that there are two consequences of increasing macroion charge: the crys- 
tallization window narrows, and it shifts to smaller values of k&I '/e. This latter effect is 
useful for protein crystallization where we might not have the freedom to raise tempera- 
ture (lest it lead to denaturation): increasing protein charge might be one way to move the 
crystallization window into the region of k^T/e accessible exprimentally. The tradeoff is 
that adding charge causes the window to narrow, but this can be counteracted by adding 
salt. Fig. 3(b) shows how, for two fixed bond strengths e/k^T, one can alter charge and salt 
concentration to widen the crystallization window. Note that this window narrows as bond 
strength increases. 

The shrinking of the crystallization window with the bond strength can be understood 
by writing the net interaction between the particles in terms of the attractive and repulsive 
components e net = — e + e rcp , where e rep = 2A/| s /z is the electrostatic repulsion energy per 
crystal contact. In terms of these variables we can express the solution-crystal coexistence 
line, given by fi = fi 3 , as 

e rep = e + ^(3A;BTln(A-l)+/io). (7) 

Similarly, we can express the solution-gel coexistence line as 

' (A; B Tln(A 3 -l)+/i ), (8) 



p A A 

where we have used A 31 to write the electrostatic free energy of the gel in terms of e rep . 

8 




(b) 






12 
8 
4 






1 1 


1 




1 




















- 


solution „--* 


e/kBT 


= 2.5 

sol 


gel 

- xtal 

-gel 


- 




gel 

1 1 


i 




1 







0.1 0.2 0.3 0.4 0.5 

c/M 



FIG. 3. (a) Phase diagram in the temperature-density plane showing solution-gel (red) and crystal- 
solution (gray) coexistence lines, for macroion charges q = (solid) and q = 6e (dotted) for salt 
concentration 0.1 M. As macroion charge increases the crystal can be rendered stable by decreasing 
temperature, but counterion entropy increasingly favors the gel over the crystal. The 'window' 
(shaded) within which the crystal is stable and the solution is stable against gelation therefore 
shrinks as charge increases, (b) Phase diagram in the charge-salt concentration plane for two bond 
strengths ejk^T . As charge increases, adding salt lessens the entropic difference between gel and 
crystal, widening the desired shaded region between the sol-gel and sol-crystal lines. Note that this 
region narrows as bond strength increases. 



Eqs. 7 and 8 are plotted in Fig. 4. We see that strongly attractive particles (large e) 
can be brought within the crystallization window by increasing e rep . However, this win- 
dow becomes narrower with increasing e, and above a critical attraction vanishes entirely. 
This happens because strongly attractive particles require large charges to weaken particle- 
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FIG. 4. Phase diagram showing the lines jjlq = n\ (dotted) and [Iq = /x 3 (solid) as a function of the 
attractive and repulsive contributions to the net interaction (here A31 = 0.6 and c = 0.04). 



particle interactions enough to destabilize the gel, but large charges in turn recruit dense 
screening layers which act to favor the gel over the crystal. This effect can be mitigated 
by working under conditions for which electrostatic interactions are approximately pairwise- 
additive (large A 31 ). Increasing the value of A 31 decreases the slope of the solution-gel 
coexistence line (dotted line in Fig. 4). This widens the crystallization window and in- 
creases the critical value of e at which the solution-gel line overtakes the solution-crystal 
line. In practice, A31 can be increased by simultaneously increasing salt concentration and 
particle charge. For example, within our model system a macroion charge q = 2.8 in 100 
mM salt, or a charge of q = 8.1 in 800 mM salt, both destabilize crystal bonds by 0.5 k-^T. 
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However, the former condition destabilizes bonds in the gel phase by only 0.17 k&T, while 
the latter condition destabilizes bonds by 0.25 k^T. Thus, high salt/high charge conditions 
are more conducive to crystallization than low salt/low charge conditions. We speculate that 
the addition of multi-valent counterions may further reduce the entropic cost of neutralizing 
the crystal, athough a quantitative analysis of this scenario is beyond the limitations of our 
PB model. 

IV. CONCLUSION 

We have studied a minimal model of the crystallization and gelation of charged attractive 
spheres in aqueous salt solution. We find that there exists an electrostatic driving force 
that favors the gel over the crystal, because there exists more ion-accessible volume for 
confined counterions around a gel than within a crystal. This effect renders the effective 
interactions between macroions non-pairwise-additive, and is responsible for a narrowing of 
the crystallization window as macroion charge is increased. 
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Appendix A: Poisson-Boltzmann free energies 

For the environmental conditions considered in this paper we find that solutions $^ of 
the Poisson-Boltzmann (PB) equation [Eq. (2) in the main text] are closely approximated 
by solutions $J] n of its linearized counterpart, the Debye-Huckel (DH) equation. We present 
numerical results derived from the PB equation, but also present, for the reader's conve- 
nience, the well-known analytic solutions of the DH equation. For completeness, we also 
consider a 2d membrane (in 2d the PB equation has an analytic solution). We draw the 
same qualitative conclusions from both levels of theory. 
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To compute the electrostatic free energy cost of assembling the building blocks into each d- 
dimensional structure, we compute the free energy per particle in the assembly, fd = F^/N 
and subtract from this the electrostatic free energy of the monomeric building blocks in 
solution. The per-particle electrostatic free energy of assembly is then Afd — fd ~ fo- To 
determine fo we solve the PB equation for the potential of an isolated sphere of radius a 
carrying charge q. The corresponding solution of the DH equation is 

Sfr(x) = ^l e —. (Al) 

1 + Ka x 

Here £b = e 2 / {^eksT) is the Bjerrum length. 

ID assembly. We treat the ID assembly as an infinitely long cylinder of radius a and linear 
charge density p = q/2a, and obtain $1 from the PB equation in plane polar coordinates. 
The corresponding solution of the DH equation is 29 

*i n (*) = ^-= 2 Q Ko(x). (A2) 

Here x is the scaled distance from the cylinder center and K n is the n th order modified Bessel 
function of the second kind. 

2D assembly. We treat each surface of the 2D assembly as an infinitely extended plane 
carrying areal charge density a = qA2/ c 2.na 2 . Here A2 — 0.91 is the area occupied by close- 
packed spheres in a sheet. The exact solution of the PB equation in planar geometry at a 
scaled distance x from the plane surface is 

. . , , /l + e- x tanh(0 o /4)\ /A x 

where 0o — 2sinh _1 [27r£BC"/( K e)] is the electrostatic potential at the plane surface. From 
Eqs. (3) and (4) in the main text we find 

nk B T A 2 . . .. 

U^A,^ — —, (A4) 



and 

8c 3A 2 -2Aln (^) 

k 1-A 



S 2 /k B = -4shcct^^ \ x2 1 A , (A5) 



where A = tanh(0 o /4) and A s h cc t is the area of the sheet. 

3D assembly. Following previous work on colloidal systems 16,19 ' 20 we treat the 3D as- 
sembly as a collection of spherical macroions of radius a, each of which is surrounded by a 
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spherical aqueous cavity of radius b. We assume that macroions are close-packed at volume 
fraction .4.3 = (a/b) 3 ~ 0.74. To calculate the electrostatic free energy of this assembly 
we assume that the electric field on the surface of each aqueous cavity vanishes 20 (valid 
for assemblies whose characteristic linear size is much greater than the Debye length). We 
also assume that the field at the macroion surface is unaffected by the presence of salt and 
counterions. We find the potential $3 from the appropriate PB equation. This solution is 
closely approximated by linearizing the PB equation around the average potential between 
macroions. By writing $3(2;) = $ 3 in (:r) + 0, where x is the scaled distance from the macroion 
center, and the mean potential 

= sinh -i ( 3 ^-—] (A6) 



\87rc a 3 (l — As) j 

is given by a jellium model 20 ' 25 , we impose the boundary conditions described above and 
find, after some algebra, 

$W ? a = ^0 x 

3 Ky) e-«+?{a + 1)03 - 1) - e a ~P{a - l)(/3 + 1) 

; + y — '- -tanhc/.. (A7) 



y V 

Here E = ge/t(cosh0) 1 / 2 / (Anek-sT); y = (cosh (ft) 1 / 2 nr; a = (cosh0) 1//2 /ta; and (3 = 
(cosh 0) 1//2 'Kb. 

Our model is likely to become unreliable at high salt concentrations and at low values 
of the packing fractions ^2,3 ■ Under such conditions the approximation of taking filament 
and sheet surfaces to be smooth becomes unrealistic because the aqueous volume of surface 
corrugations becomes comparable to the total volume of the screening layer (accounting for 
steric corrections to the PB equation 30 then becomes necessary). In order to neglect these 
corrugations, it is necessary that both the Debye length and the Gouy-Chapman length, 
^gc — e/27r£B<7, exceed the characteristic length scale of the surface cavities. For this reason 
we restrict our analysis below to salt concentrations below c = 0.5 M and q = 10, for our 
chosen macroion radius a = 1.6 nm. 
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